What if we wanted to look at population in 20 years given an initial condition
Two options
explicit solution to differential equation is known; e.g. you can integrate both sides of the equation! Not always possible but lets look at a case where it is possible
must be solved by iteration; this is what we do when we can’t integrate both sides
## function (T, P0, r, K)
## {
## P = P0 * exp(r * T)
## if (P > K) {
## P = K
## }
## return(P)
## }
# gives population after any time given an initial population
# 20 rabbits, growth rate of 0.01 how many in 30 years
exppop(T=30, P0=20, r=0.01, K=1000)## [1] 26.99718
# if we want to see how population evolves over time - generate a time series by running our model for each point in time
initialrabbits = 20
years = seq(from=1, to=100, by=2)
Ptime = years %>% map_dbl(~exppop( P0=initialrabbits, r=0.01, K=1000, T=.x))
# keep track of what times we ran
Ptime = data.frame(P=Ptime, years=years)
ggplot(Ptime, aes(years,P))+geom_point()+labs(x="years",y="Rabbit Population")Population models can be discrete (rather than continuous)
So we could implement them as difference equations and iterate
source("../R/discrete_logistic_pop.R")
# notice how a for loop is used to iterate
# how many rabbits after 50 years given a growth of 0.1
# starting with 1 rabbit - but a carrying capcity of 500
discrete_logistic_pop## function (P0, r, K, T = 10)
## {
## pop = P0
## for (i in 1:T) {
## pop = pop + r * pop
## pop = ifelse(pop > K, K, pop)
## }
## return(pop)
## }
## [1] 11.4674
## [1] 12.18249
## [1] 12.18249
## [1] 11.4674
# why are they different
# look at trajectories
growth_result = data.frame(time=seq(from=1,to=100))
growth_result$Panalytic = growth_result$time %>% map_dbl(~exppop( P0=1,r=0.05, K=200,T=.x ))
growth_result$Pdiscrete = growth_result$time %>% map_dbl(~discrete_logistic_pop( P0=1,r=0.05, K=200,T=.x ))
tmp = growth_result %>% gather(key="Ptype",value="P",-time)
ggplot(tmp, aes(time,P, col=Ptype))+geom_point()Using a solver….when you can’t do the integration by hand :)
For example, if you added a carrying capacity as threshold where it stops growing
In this case we integrate by iteration - approximates analytic integration
Implement the differential equation as a function that
*inputs time, the variable(s) and a parameter list
(it needs time even though you don’t use it)
My convention: name derivative functions starting with d to remind myself that they are computing a derivative
Only works for Ordinary Differential Equations - single independent variable (in our case time)
Partial differential equations - more than 1 independent variable (e.g x and y if changing in space)
R has a solver called ODE for solving ordinary differential equations from package desolve
## function (time, P, r)
## {
## dexpop = r * P
## return(list(dexpop))
## }
library(deSolve)
initialrabbits = 20
years = seq(from=1, to=100, by=2)
# run the solver
Ptime = ode(y=initialrabbits, times=years, func=dexppop,parms=c(0.01))
head(Ptime)## time 1
## [1,] 1 20.00000
## [2,] 3 20.40404
## [3,] 5 20.81623
## [4,] 7 21.23675
## [5,] 9 21.66576
## [6,] 11 22.10344
colnames(Ptime)=c("year","P")
# notice that there are additional pieces of information year, including the method used for integration
attributes(Ptime)## $dim
## [1] 50 2
##
## $dimnames
## $dimnames[[1]]
## NULL
##
## $dimnames[[2]]
## [1] "year" "P"
##
##
## $istate
## [1] 2 52 105 NA 6 6 0 36 21 NA NA NA NA 0 1 1 NA NA NA
## [20] NA NA
##
## $rstate
## [1] 2.00000 2.00000 99.98839 0.00000 1.00000
##
## $lengthvar
## [1] 1
##
## $class
## [1] "deSolve" "matrix"
##
## $type
## [1] "lsoda"
# this also means you need to extract just the data frame for plotting
ggplot(as.data.frame(Ptime), aes(year,P))+geom_point()+labs(y="Population", "years")# this also works (of course function can be by order)
Ptime=ode(initialrabbits, years, dexppop,0.03)
colnames(Ptime)=c("year","P")
ggplot(as.data.frame(Ptime), aes(year,P))+geom_point()+labs(y="Population", "years")You can play a bit with changing your function to something that you can’t integrate “by hand”
BUT we might want more parameters
to work with ODE, parameters must all be input as a single list; similar to how we return multiple outputs from a function
see example below..lets add a carrying capacity
## function (time, P, parms)
## {
## dexpop = parms$r * P
## dexpop = ifelse(P > parms$carry_capacity, 0, dexpop)
## return(list(dexpop))
## }
# create parameter list
initalrabbits=2
newparms = list(r=0.03, carry_capacity=300)
#apply solver
results=ode(initialrabbits, years, dexppop_play,newparms)
head(results)## time 1
## [1,] 1 20.00000
## [2,] 3 21.23675
## [3,] 5 22.54993
## [4,] 7 23.94435
## [5,] 9 25.42499
## [6,] 11 26.99719
# add more meaningful names
colnames(results)=c("year","P")
# plot
ggplot(as.data.frame(results), aes(year,P))+geom_point()+labs(y="Population", "years")# try again with different parameters
alternativeparms = list(r=0.04, carry_capacity=500)
results2=ode(initialrabbits, years, dexppop_play,alternativeparms)## DLSODA- At current T (=R1), MXSTEP (=I1) steps
## taken on this call before reaching TOUT
## In above message, I1 = 5000
##
## In above message, R1 = 81.5108
##
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps
## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
## time 1
## [1,] 1 20.00000
## [2,] 3 21.66575
## [3,] 5 23.47022
## [4,] 7 25.42499
## [5,] 9 27.54257
## [6,] 11 29.83650
colnames(results2)=c("year","P_parm2")
# plot
ggplot(as.data.frame(results2), aes(year,P_parm2))+geom_point()+labs(y="Population", "years")# compare by combining into a single data frame
both = inner_join(as.data.frame(results), as.data.frame(results2))## Joining with `by = join_by(year)`
both_p = both %>% gather(key="model", value="P", -year)
ggplot(both_p, aes(year,P, col=model))+geom_point()+labs(y="Population", "years")Remember we have 3 ways now to calculate population
analytical solution - based on integration (exppop.R) BEST
using an ode solver for numerical approximation (exppop_play.R)
numerical integration using in discrete steps (discrete_logistic_pop.R)
Finally look at continuous derivative using ODE solve Needs initial conditions differential equation *parameters
## function (time, P, parms)
## {
## dexpop = parms$r * P
## dexpop = ifelse(P > parms$carry_capacity, 0, dexpop)
## return(list(dexpop))
## }
# set up using the same parameters
pcompare = list(r=r,carry_capacity=K)
# now run our ODE solver
result = ode(P0, growth_result$time, dexppop_play, pcompare)
head(result)## time 1
## [1,] 1 1.000000
## [2,] 2 1.051273
## [3,] 3 1.105172
## [4,] 4 1.161835
## [5,] 5 1.221404
## [6,] 6 1.284027
# we already have time - so just extract population
growth_result$Pdifferential = result[,2]
# compare all 3 approaches
tmp = growth_result %>% pivot_longer(cols=-time,names_to="Ptype",values_to="P")
ggplot(tmp, aes(time,P, col=Ptype))+geom_point()All differential and difference equations are approximations The analytical solution is exact
Notice that differential equations is a bit more accurate!
Diffusion can be implemented as a partial differential equation
Complicated to solve - but there are tool in R for specific types of partial differential equations
More info on differential equations in R
Diffusionn would require partial derivatives - time and space! it gets much more tricky …beyond this course
But we can appoximate diffusion with a difference equation - and iterative to get an estimate of how diffusion works
Example of Diffusion - difference equation implementation to see what some issues can be
source("../R/diffusion.R")
# run our diffusion model (iterative difference equation) with initial concentration of 10, for 8 timestep (size 1m), and 10 space steps (size 1s)
# using diffusion parameters 0.5 s/m2, 10 m2
result = diff1(initialC=10, nx=10, dx=1, nt=8, dt=1, D=0.5, area=10)
# a list is returned with our 3 data frames for concentration (conc), qin and qout
result## $conc
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 10.000000 0.000000 0.000000 0.0000000 0.0000000 0.000000000 0.000000000
## [2,] 7.500000 2.500000 0.000000 0.0000000 0.0000000 0.000000000 0.000000000
## [3,] 6.250000 3.125000 0.625000 0.0000000 0.0000000 0.000000000 0.000000000
## [4,] 5.468750 3.281250 1.093750 0.1562500 0.0000000 0.000000000 0.000000000
## [5,] 4.921875 3.281250 1.406250 0.3515625 0.0390625 0.000000000 0.000000000
## [6,] 4.511719 3.222656 1.611328 0.5371094 0.1074219 0.009765625 0.000000000
## [7,] 4.189453 3.142090 1.745605 0.6982422 0.1904297 0.031738281 0.002441406
## [8,] 3.927612 3.054810 1.832886 0.8331299 0.2777100 0.064086914 0.009155273
## [,8] [,9] [,10]
## [1,] 0.0000000000 0 0
## [2,] 0.0000000000 0 0
## [3,] 0.0000000000 0 0
## [4,] 0.0000000000 0 0
## [5,] 0.0000000000 0 0
## [6,] 0.0000000000 0 0
## [7,] 0.0000000000 0 0
## [8,] 0.0006103516 0 0
##
## $qout
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 25.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000 0.000000000
## [2,] 12.500000 6.250000 0.000000 0.000000 0.00000000 0.00000000 0.000000000
## [3,] 7.812500 6.250000 1.562500 0.000000 0.00000000 0.00000000 0.000000000
## [4,] 5.468750 5.468750 2.343750 0.390625 0.00000000 0.00000000 0.000000000
## [5,] 4.101562 4.687500 2.636719 0.781250 0.09765625 0.00000000 0.000000000
## [6,] 3.222656 4.028320 2.685547 1.074219 0.24414062 0.02441406 0.000000000
## [7,] 2.618408 3.491211 2.618408 1.269531 0.39672852 0.07324219 0.006103516
## [8,] 0.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000 0.000000000
## [,8] [,9] [,10]
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 0 0 0
## [4,] 0 0 0
## [5,] 0 0 0
## [6,] 0 0 0
## [7,] 0 0 0
## [8,] 0 0 0
##
## $qin
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0 25.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000
## [2,] 0 12.500000 6.250000 0.000000 0.000000 0.00000000 0.00000000
## [3,] 0 7.812500 6.250000 1.562500 0.000000 0.00000000 0.00000000
## [4,] 0 5.468750 5.468750 2.343750 0.390625 0.00000000 0.00000000
## [5,] 0 4.101562 4.687500 2.636719 0.781250 0.09765625 0.00000000
## [6,] 0 3.222656 4.028320 2.685547 1.074219 0.24414062 0.02441406
## [7,] 0 2.618408 3.491211 2.618408 1.269531 0.39672852 0.07324219
## [8,] 0 0.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000
## [,8] [,9] [,10]
## [1,] 0.000000000 0 0
## [2,] 0.000000000 0 0
## [3,] 0.000000000 0 0
## [4,] 0.000000000 0 0
## [5,] 0.000000000 0 0
## [6,] 0.000000000 0 0
## [7,] 0.006103516 0 0
## [8,] 0.000000000 0 0
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 10.000000 0.000000 0.000000 0.0000000 0.0000000 0.000000000 0 0 0
## [2,] 7.500000 2.500000 0.000000 0.0000000 0.0000000 0.000000000 0 0 0
## [3,] 6.250000 3.125000 0.625000 0.0000000 0.0000000 0.000000000 0 0 0
## [4,] 5.468750 3.281250 1.093750 0.1562500 0.0000000 0.000000000 0 0 0
## [5,] 4.921875 3.281250 1.406250 0.3515625 0.0390625 0.000000000 0 0 0
## [6,] 4.511719 3.222656 1.611328 0.5371094 0.1074219 0.009765625 0 0 0
## [,10]
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 0
## [5,] 0
## [6,] 0
# or if you prefer this orientation (Distance on x axis)
filled.contour(t(result$conc), ylab="Time", xlab="Distance")# changes diffusivity and other parameters particularly
# diffusivity, dx and dt
res=diff1(initialC=10,nx=10,dx=1,nt=10,dt=30,D=0.006,area=1)
filled.contour(res$conc, xlab="Time", ylab="Distance")# we can also see how much material moved from place to place each time step
filled.contour(res$qin, xlab="Time", ylab="Distance")Try running the diffusion model with different time steps, space steps and parameters
# what if we increase diffusivity
resfast=diff1(initialC=10,nx=10,dx=0.5,nt=10,dt=10,D=0.08,area=1)
filled.contour(resfast$conc, xlab="Time", ylab="Distance")# Discretization Issue Example
resunstable=diff1(initialC=10,nx=10,dx=1,nt=10,dt=10,D=0.8,area=1)
filled.contour(resunstable$conc, xlab="Time (fraction of hour)",ylab="Distance Along Path (m)", main="Pollutant Diffusion")# this illustrates the problem with difference equations (and the challenges that methods for numerical integration try to overcome)
# if things are changing quickly we need to use much smaller time, space steps to avoid overshoot and instability
# so lets cut our step size by 10 (dt) (but then add 10 more steps (nx) to cover the same distance)
resunstable=diff1(initialC=10,nx=100,dx=1,nt=10,dt=1,D=0.8,area=1)
filled.contour(resunstable$conc, xlab="time",ylab="Distance Along Path")